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Preface 

The  procedure  for  an  exact  test  of  independence  in  general  m  x  n 
contingency  tables  has  been  available  for  years.  However,  as  with 
much  of  our  knowledge,  the  application  of  theory  in  a  practical  sense 
has  had  to  wait  for  the  developement  of  high  speed  computers.  The 
program  presented  here  should  be  most  useful  in  those  instances 
where  a  test  by  approximation  (although  relatively  simple  and  easily 
calculated)  is  either  questionable  or  clearly  inadvisable. 

I  wish  tc  ‘■hark  Dr.  David  R.  Barr,  my  thesis  advisor,  for  sug¬ 
gesting  the  topic  and  for  his  infinite  patience  in  prodding  me  to  the 
final  solution. 

Marvin  F.  Schwartz,  Jr. 
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Abstract 

General  m  x  n  contingency  tables  are  discussed  along  with  re¬ 
lationships  common  to  all  such  tables.  The  background  for  various 
probability  models  is  presented  as  well  as  the  theory  for  a  general 
conditional  test.  The  general  test  is  then  applied  specifically  to  the 
probability  model  selected  followed  by  a  numerical  example. 

A  brief  history  of  past  efforts  in  computing  the  exact  test  is 
provided.  Finally,  the  major  portions  of  the  FORTRAN  program 
are  presented  along  with  some  illustrative  examples.  The  final 
program  is  included  as  a  useful  entity. 


Results  of  Results  of 

1st  Ejection  2nd  Ejection 


Group 

Number 

Non-injured 

Injured 

Fatal 

Non-injured 

84 

43 

27 

14 

Injured 

35 

22 

12 

1 

Totals 

119 

65 

'  39 

15 

(From  Ref  28) 

This  data  represents  the  total  United  States  Air  Force  experience 
with  multiple  ejections  and  as  arranged  forms  the  well  known  con¬ 
tingency  table.  Reference  will  be  made  back  to  this  data  later  to 
show  the  analysis  by  approximation  as  reported  and  then  the  analysis 
by  an  exact  test. 

General  M  X  N  Contingency  Tables 


A  FORTRAN  PROGRAM  FOR  AN 
EXACT  TEST  OF  INDEPENDENCE 
IN  AN  M  X  N  CONTINGENCY  TABLE 

I.  Introduction 


The  May,  1970  issue  of  Aerospace  Medicine  contains  an  article 


which  reports  the  following  data  on  pilots  who  have  made  multiple 


ejections: 


Table  I 

USAF  Experience  With  Multiple  Ejections 
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Suppose  that  K  individuals  arc  classified  according  to  two 
criteria  X  and  V.  Suppose  further  that  there  are  m  classifications 
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for  X  and  n  classifications  for  Y.  A  single  individual  would  then 
be  classified  according  to  Xj  (i  =  1,2,...,  m)  and  Y.  (j  =  1, 2,  . .  .  ,  n). 
If  the  number  of  individuals  belonging  to  both  Xj  and  Yj  is  Ajj  we 
have  an  m  x  n  contingency  table.  Such  a  table  is  shown  in  Figure  1. 


Y1 

Y2 

•  •  • 

•  •  • 

Y 

n 

Totals 

X1 

A 

11 

A 

12 

D 

B 

Aln 

R1 

X2 

A21 

A22 

m 

A2n 

R2 

• 

■ 

• 
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• 

• 

• 

• 

■ 

• 

El 

B 

• 

• 

• 

• 

• 

• 

• 

• 

s 

X 

m 

Aml 

Am2 

B 

B 

Amn 

R 

m 

Totals 

C2 

B 

B 

c 

n 

K 

Figure  1.  General  m  x  n  Contingency  Table 
Some  Relationship  s 

There  are  several  relationships  common  to  all  contingency 

tables  and  essential  to  a  discussion  of  the  subject.  If  row  totals 

are  designated  as  R.  (i  =  1,2 . m)  then  R.  =  )  "  .  A...  Sim- 

i  1  Lj=l  ij 

ilarly,  if  column  totals  are  designated  as  Cj  (j  =  1, 2, . .  . ,  n)  then 

_  V-1  m  t~»m  n 

C .=  )  .  A...  Further,  )  R.  =  K,  )  .  C.  =  K,  and 

J  L-,  1=1  lj  l—t  l-l  l  La}~l  J 

Em  nn 

i=l/  t  j=i  A.j  =  K,  where  K  is  the  total  sample  size.  Let  p..  be 
the  probability  that  an  event  falls  in  the  cell  corresponding  to  the 
ith  row  and  jth  column.  Let  p.  be  the  probability  that  an  event 
falls  into  the  ith  row  and  q.  be  the  probability  that  an  event  falls 

J 
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into  the  jth  column.  Then  Pj=  ^  j*-i  Pjj'  *^j=  ^  i~l  Pij  anc* 

Em  on  c— >m  r— i  n 

._.)  p..  =  1.  Also,  )  .  ,  p.  =  1  and  )  q.  =  1. 

i=l£_,j=l  Hij  Z_.  i=  1  *i  Li  J=1  J 


Extension  to  Higher  Dimensions 

Although  the  material  presented  here  deals  with  just  two  criteria 
the  only  limit  to  the  number  of  criteria  which  can  be  dealt  with  is 
one  of  practicality.  If,  for  example,  the  elements  of  a  population 
were  to  be  classified  according  to  three  criteria  X,  Y,  and  Z  with 

classifications  Xi  (i  =  1,  2,  ....  r),  Yj  (j  =  1,2 . s),  and 

(k  =  1,2 . t),  a  sample  of  K  individuals  could  be  classified  in 

a  three-way  r  x  s  x  t  contingency  table.  Thus  there  is  no  theo¬ 
retical  limit  to  the  number  of  possible  criteria  and  the  treatment 
of  the  theory  with  added  dimensions  poses  no  special  problems. 
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II.  Testing  For  Independence 

The  Hypothesis  of  Independence 

The  contingency  table  is  merely  a  convenient  way  to  arrange 
data.  The  data  itself  can  lead  to  numerous  questions  of  cause  and 
effect,  or  interaction  (Ref  23:260).  The  question  to  be  addressed 
here  is  one  of  no  interaction,  or  independence.  That  is,  is  class¬ 
ification  by  one  criterion  independent  of  classification  by  the  other 
criterion?  Referring  back  to  the  general  m  x  n  table  of  Figure  1 
we  can  ask  if  classification  in  any  column  is  independent  of  the 
row  entry'  of  an  individual.  If  the  entries  are  independent  then  the 
probability  of  an  individual  being  in  a  particular  cell  is  merely  the 
product  of  the  probabilities  of  the  row  and  column  associated  with 
that  cell.  The  null  hypothesis  can  therefore  be  formulated  as 

Ho:  pij  =  piqj  l'  z  1’2 . .  j  =  1,2,  ...  ,n). 

The  probability  model  which  will  be  used  here  is  that  which 
G.  A.  Barnard  named  as  '  j  Independence  Trial  (Ref  2)  which 
treats  the  table  under  consideration  as  though  the  marginal  totals 
(i.e.  row  and  column  totals)  are  fixed  in  advance.  Other  models 
are  the  Double  Dichotomy  which  treats  the  marginal  totals  as  ran¬ 
dom  variables  and  the  Contingency  Trial  (or  Homogeneity)  which 
considers  one  set  of  marginals  as  fixed  and  the  other  set  as  random 
variables  (Ref  2,  3,  4,  17:549-554;  27,  29,  33).  The  Independence 
Trial  model  is  selected  here  since  it  is  the  only  model  which  pro- 
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vides  for  an  exact  test  and,  from  among  the  three  models  avail¬ 
able,  it  provides  the  most  powerful  test  (Ref  17:554;  29).  Thus, 
by  fixing  the  marginal  totals,  we  will  be  using  a  conditional  test 
based  on  the  probability  of  a  given  table  having  fixed  marginal 
totals.  As  will  be  shown,  this  conditional  distribution  is  known 
and  permits  proceeding  with  an  exact  test. 

In  the  vernacular  of  statistical  hypothesis  testing  what  is  de¬ 
sired  is  a  "two-tailed"  test.  Had  the  hypothesis  been  one  which 
asked  if  entry  in  a  particular  cell  was  more  or  less  likely  to  occur 
as  a  result  of  the  associated  row  and  column  we  would  then  have 
sought  a  "one-tailed"  test.  Here  we  are  trying  to  determine 
(within  limits)  whether  classification  by  one  criterion  is  unaffected 
in  either  direction  (high  or  low)  by  the  other  criteria  (Ref  23:172). 
Hence,  we  require  a  "two-tailed"  or  non-directional  test. 

The  Background  for  a_  Conditional  Test 

To  construct  a  conditional  test,  consider  the  general  situation 

of  a  random  vector  X  =  (X.,  X  , . .  .  ,  X  )  with  a  distribution  f(X;e) 

■*  1  2  n  — 

and  involving  the  unknown  parameters  e  z  (e  j ,  •  •  •  >  es)’  ^  $ 

A 

has  a  set  of  joint  sufficient  statistics  0,  then  the  joint  density  of 

A 

X  and  0  may  be  written: 

f(x,5;e)=  fi(x|«)f2(S;e)  (1) 

*  A 

Now  the  conditional  distribution  of  X  given  ©  =  e  does  not  depend 
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on  e  because  the  ©  are  joint  sufficient  statistics  and  for  each  value 
of  5  we  can  find  numbers  A(e-)  and  B(6j),  (i  =  1,  2, ....  s)  which 


are  functions  of  e.,  so  that 


Ja,Si,  •  •  •  Ja,s,  m 

i  s 

where  (X  is  some  desired  significance  level  for  a  critical  region 
C.  Then  with  some  sample  space  S  having  ar.  acceptance  region 
A  we  have 

r  .  -1  r°°  (-00  r  B(6j)  rB(es) 

p[a(©)«X«:  B(@)lr;J  ...  /  . ../  .  f .  (xle)f2(e;e)dxd©  (3) 

L  -  ~  -Uoo  J-ooJ  A(V  J  M*s)  - 

rOO 

=  [  •••/  H*aK2L5Lsld! 

-Loo  J-OO 

rOO  roo 

=d  -a)/  ••■/ 

-Loo 


=  i  -a 


When  X  is  a  single  variable  and  9  a  single  statistic  the  acceptance 
region  A  is  as  shown  in  Figure  2. 
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Thus  there  exists,  for  each  value  of  e,  an  interval  (A(e),B(e))  such 
that  P^A(e)«  X«B(e)|©  =  &J  =  1  -0£.  Naturally,  when  this  concept 
is  extended  to  the  multi-parameter  case  a  geometric  represen¬ 
tation  becomes  impossible  (Ref  25:258;  30:525). 

Notice  that  this  developement  is  dependent  on  formula  (1)  and 
requires  the  joint  sufficient  statistics  ©  in  order  to  make  the  cal¬ 
culations  shown.  If  f  j(x|8)  were  not  free  of  e,  then  calculation 
of  the  probability  would  depend  on  unknown  parameters  and  the 
procedure  used  above  would  be  inappropriate. 


The  Test  as  Employed 

The  hypothesis  of  independence  was  stated  as,  Hq:  p.^  r  P^j* 

Unfortunately,  the  parameters  p.  and  are  unknown.  However, 

it  can  be  shown  that  the  estimators  R^/K  (i  =1,2,...,  m)  and 

C./K  (j  =  1, 2, . . . ,  n)  form  a  set  of  joint  sufficient  statistics  for 
J 

p^  and  q^  respectively.  By  a  one  for  one  transformation,  Rj 
(i  =  1, 2, ....  m)  and  (j  =  1,  2, . . . ,  n)  are  also  sufficient  statis¬ 
tics  (Ref  25:171)  and  we  can  determine  the  exact  probability  for 
any  contingency  table  by  means  of  the  following  formula: 

,k 


PL  = 


-I— Tk  -i— rr 

TT  a™i  ! 

|m=l|  1  i  =1  'm 


(4) 


where  is  the  conditional  probability  of  an  (rj  x  r^  x*..x  r^) 
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K-variate  contingency  table  L  =  (1.  .  .  ),  where  1;  ;  :  i 

8  V2’  k  12'  ‘  k 

the  entry  in  the  cell  identified  by  the  K  indices  i  given  all  the 
marginal  totals 

_  yTi  yrm-i  v"rm+i  yrk 

amim"  ^il=1‘“Zjim-l=1^im+i::1“,^ik=1 

This  formula  was  first  presented  by  Freeman  and  Halton  (Ref  10) 

and  a  more  recent  derivation  is  presented  by  Halton  (Ref  13). 

Eq  (4)  appears  in  more  restricted  forms  in  the  literature  and  for 

use  here  with  just  two  criteria  can  be  reduced  to: 


(5) 


TTuffc,  r>!cJ! 

~  kinr=.n;=1  v 

using  the  notation  first  presented.  Numerous  derivations  are  avail¬ 
able  for  Eq  (6).  The  derivation  for  a  2  x  2  table  was  first  published 
by  Fisher  in  1934  (Ref  7).  Irwin  derived  it  independently  in  1935 
and  the  test  as  applied  to  a  2  x  2  table  is  often  referred  to  as  "The 
Fisher -Irwin  Exact  Test".  Most  derivations  are  based  on  a  like¬ 
lihood  ratio  test  statistic  (Ref  25:313-316).  The  right  side  of  Eq  (6) 
is  the  multivariate  hyper  geometric  distribution  and  gives  the  condi¬ 
tional  distribution  of  the  A.,  given  R.  and  C;.  Since  the  last  cell 
entry  in  each  row  can  be  determined  from  the  first  (n-1)  entries  in 


PL  = 


(6) 


that  row  and  the  row  total,  and  the  last  cell  entry  in  each  column  is 
fixed  by  the  first  (m-1)  entries  and  the  column  total,  the  distri¬ 
bution  involves  (m-1)  x  (n-1)  independent  variables.  Eq  (6)  can 
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thus  be  substituted  into  a  discrete  form  of  Eq  (2)  to  determine  an 
appropriate  acceptance  region.  In  applying  the  exact  test  to  a 
general  m  x  n  table  Freeman  and  Halton  (Ref  10)  have  suggested 
the  following  rule  for  computing  an  exact  P  value  (Ref  12).  The 
conditional  probability  of  an  observed  table  can  be  computed  from 
Eq  (6).  The  conditional  probability  of  observing  any  other  table 
having  the  same  marginal  totals  is 

__m  n 

P(x=„|R=r,  ^..TWW.'V 

.  ‘TTS,TT;,v 

where  y™  y  ?  x..  =  k.  y™.  =  c.,  andyn_.  x..  r  r..  We 

Z-.i=lZ_iJ=l  ij  Z_.1=1  J  L>  J=1  >J  » 

consider  all  possible  configurations  of  such  x's  such  that 

P(X=x|R=r,  Crc)  £  P(A=a|R=r,  Crc)  ( 

Then  P  =  J]  P(X.rx.|Rrr,  Crc)  where  the  sum  is  with  respect  to 
it I  -i  l  ~  -  ~  ~ 

the  x's  subject  to  the  constraint  (8). 

Kendall  and  Stuart  (Ref  17:552)  give  a  numerical  illustration 
using  the  following  table: 


Table  II 
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Given  fixed  marginal  totals  the  probability  that  A^  s  4  is 

cl!c2!rl!r2!  5 !  37 !20 ! 22  ! 

k!all!al2!a21!a22!  =  4Z!4!16!1  !21!  = 

Probabilities  for  all  possible  X  given  r,r,c,c  are  plotted 

11  1  £*  1  Ct 

in  Figure  3. 


Figure  3.  P(Xjj=  Xjj)  versus  Xjj 


P  -  4_»  P(X, ,  =  x, ,  |R=r,  C=c)  where  the  sum  is  over  the  x's 
iCI  llj  11. - - 

subject  to  the  constraint  {8)  is  .  1745.  Thus  at  a  significance  level 
of  .  05  or  .  10,  for  example,  we  would  fail  to  reject  the  hypothesis 
of  independence  for  the  given  table.  The  distribution  is  a  univariate 
hypergeometr ic  and  is  unimodal  while  the  test  is  non-directional. 

An  observed  table  with  the  same  marginal  totals  but  with  Ajj  equal 
to  0  or  5  would  have  led  to  rejection  of  the  hypothesis  at  a  signif¬ 
icance  level  of  .05. 
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III.  The  Need  for  a  Computer  Program 
Chi-aquare  and  its  Limitationa 

Aa  can  be  aeen,  the  formula  for  computing  an  individual  prob¬ 
ability  can  lead  to  an  extremely  laborioua  calculation.  In  worda. 


_  Product  of  all  (border  totala!) 
L”  (K!)  x  product  of  all  (cell  totala!) 


(9) 


In  all  but  the  moat  trivial  problems  the  time  required  to  perform 
the  calculation  without  resorting  to  computers  could  quickly  be¬ 
come  prohibitive.  The  common  recourse  is  to  fall  back  on  an 
approximation  of  one  form  or  another.  For  example,  it  can  be 
shown  that  for  large  samples  if  the  distribution  (6)  is  replaced  by 
its  multivariate  normal  approximation,  the  criterion 


K-iK,vk,] 

"=L  -  - 


R.C./k 
*  J 


(10) 


approaches  the  chi-square  distribution  asympotically  (Ref  25). 
Wilks  (Ref  32)  derives  a  similar  statistic  with  an  approximate 
chi-square  distribution.  However,  the  chi-square  approximation 
is  limited  in  its  applicability  and  various  rules  of  thumb  prevail 
concerning  its  use  with  most  based  on  the  expected  frequency  for 


any  A.^  as  computed  by 


E  =  R.C./k 
»  J 


(ID 
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Moat  commonly,  the  approximation  is  discouraged  when  the  expected 
frequency  in  any  cell  is  less  than  five  (Ref  11:295).  Kenney  and 
Keeping  (Ref  18:307)  recommend  the  exact  method  when  the  smallest 
expected  frequency  in  any  cell  is  less  than  10  while  Mood  (Ref  25:413) 
considers  the  approximation  "fairly  good"  in  a  2  x  2  table  if  both 
column  totals  exceed  10.  Pierce  (Ref  26:115)  remarks  that,  ". .  the 
computational  convenience  of  the  chi-square  model  is  no  longer  ?. 
plausible  excuse  for  the  2x2  case. . .  " 

Other  Limited  Aids 

Exact  treatment  of  the  2x2  case  has  certainly  received  a  great 
deal  of  attention.  Best  known,  perhaps,  is  Lieberman  and  Owen's 
Tables  of  the  Hyper  geometric  Probability  Distribution  (Ref  21). 

These  were  preceded  by  Mainland  and  Murray  (Ref  24)  and  Finney 
et^al  (Ref  6)  have  provided  a  more  extensive  set.  Still  another  set 
has  been  constructed  by  Clark  (Raf  4).  Leslie  (Ref  19)  has  written 
up  a  procedure  using  binomial  coefficients  which  is  valid  for  both 
one  and  two-tailed  tests. 

Beyond  the  very  limited  case  of  a  2  x  2  table  much  less  has 
been  done.  Leyton  (Ref  20)  treats  the  2x3  table  using  Leslie's 
approach  of  binomial  coefficients  and  the  critical  region  suggested 
by  Freeman  and  Halton  (Ref  10).  Pierce  (Ref  26)  discusses  exten¬ 
sion  beyond  the  2x2  case  and  has  written  an  ALGOL  and  a  BASIC 
language  computer  program  for  the  2x3  case.  As  regards  the 
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BASIC  language  program  the  prime  algorithm  is  correct  and  the 
point  probability  computed  for  the  observed  table  is  accurate. 
However,  the  cumulative  probabilities  are  inaccurate  and  should 
not  be  relied  on.  The  program  finds  and  stores  exact  probabil¬ 
ities  for  various  possible  arrays.  The  stored  probabilities  are 
then  scanned  in  the  order  in  which  they  were  computed  until  the 
first  time  a  probability  is  found  which  exceeds  the  point  probabil¬ 
ity  of  the  observed  table.  The  program  is  then  immediately  termi¬ 
nated  with  an  output  which  is  purported  to  be,  "The  probability 
of  a  distribution  i.'s  probable  or  less  probable  than  that  observed." 
Since  the  program  terminates  abruptly  and  does  not  consider  all 
of  the  probabilities  in  the  conditional  sample  space  some  probabil¬ 
ities  ". . .  less  probable  than  that  observed. . . "  are  usually  lost 
and  not  reported.  No  theory  or  explanation  is  given  for  doing  this 
and  yet  the  program  is  classed  as  "non -directional". 

The  program  reported  here  looks  at  all  probabilities  according 
to  Eq  (6)  and  is  applicable  to  any  size  m  x  n  table  limited  only 
by  computer  capabilities. 

A  Questionable  Finding 

Consider  again  the  data  on  multiple  ejections  cited  earlier. 

For  convenience  the  data  is  repeated  below.  The  researcher  report¬ 
ed  that,  "However,  out  of  84  individuals  who  had  received  no  report- 
able  injury  on  the  first  ejection,  there  were  14  fatalities  on  the 
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Table  III 

USAF  Experience  Y/ith  Multiple  Ejections 


Results  of  Results  of 

1st  Ejection  2nd  Ejection 


Group 

Number 

Non-injured 

Injured 

Fatal 

Non-injured 

84 

43 

27 

14 

Injured 

35 

22 

12 

1 

Total 

119 

65 

39 

15 

(From  Ref  28) 

second  ejection.  A  chi-square  test  applied  to  this  finding  reveal¬ 
ed  it  to  be  highly  significant  at  the  95%  level.  " 

If  one  were  to  test  for  independence  of  classification  for  this 
table  Eq  (11)  could  be  applied  to  obtain  the  expected  frequency  in 
the  A^3  cell.  Here,  E  =  (35)  x  ( 1 5) / 1 19  :4.4  and  a  chi-square 
approximation  is  questionable. 
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IV.  Developing  the  Program 

Basic  Considerations 

In  developing  the  program  the  basic  considerations  were: 

1.  The  probability  of  the  given  array  must  be  established  and 
all  other  probabilities  in  the  sample  space  compared  with  it. 

2.  Each  feasible  array,  given  the  fixed  marginal  totals,  must 
be  determined. 

3.  The  probabilities  associated  with  each  feasible  array  must 
be  computed  and  compared  with  the  probability  of  the  given  array. 

Each  of  these  considerations  is  discussed  in  detail  and  summa¬ 
rized  by  a  numerical  example  to  clarify  the  discussion. 

The  Probability  of  the  Given  Array 

The  computation  for  the  probability  of  the  given  array  should 

employ  Eq  (6).  Since  the  result  is  a  probability  we  know  that 

0  *  P  If  1.  Although  the  equation  appears  as  a  simple  fraction, 

L 

the  number  of  factors  in  the  numerator  is  always  less  than  the  num¬ 
ber  of  factors  in  the  denominator,  each  factor  is  a  factorial,  and, 
provision  must  be  made  for  the  fact  that  0!  —  1!  =  1.  An  obvious 
approach  to  the  computation  is  to  determine  the  individual  values 
for  both  the  numerator  and  denominator  and  find  their  quotient. 

This  procedure  could  quickly  overflow  in  many  computers.  Con¬ 
sider  just  the  factor  k!  in  the  denominator  of  the  expression. 

18 

20!  =  2.43  x  10  and  it  is  obvious  that  some  extraordinarily  large 
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values  could  appear  long  before  the  final  division  of  numerator  by 
denominator.  Since  the  sum  of  all  of  the  marginal  totals  (appearing 
in  numerator)  equals  the  sum  of  each  cell  entry  plus  the  sample 
size  a  scheme  was  sought  which  would  move  from  numerator  to 
denominator  in  such  a  way  as  to  keep  any  intermediate  result  as 
small  as  possible.  The  scheme  decided  upon  requires  that  the 
array  be  arranged  as  follows: 

1.  If  the  array  is  not  square  the  number  of  rows  must  be  less 
than  the  number  of  columns. 

2.  From  top  to  bottom,  each  row  total  must  be  less  than  or 
equal  to  the  next  succeeding  row  total. 

3.  From  left  to  right,  each  column  total  must  be  less  than  or 
equal  to  the  next  succeeding  column  total. 

The  algorithm  takes  the  final  form  of: 

P  =  P<£><!>  02) 


starting  with  P  —  1. 

a  represents  the  row  totals.  It  is  set  at  one,  incremented  by 
one  up  to  the  value  of  r  j,  reset  to  one  and  incremented  up  to  the 
value  of  r^,  etc. 

b  represents  the  column  totals.  It  is  set  to  one  and  increased 
by  one  after  each  n  repetitions  (where  n  r  number  of  columns)  un¬ 
til  achieving  a  value  of  Cj.  It  is  then  increased  after  n-1  repcti- 
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tions  with  the  number  of  repetitions  between  increments  being  re¬ 
duced  each  time  b  takes  on  the  value  of  the  succeeding  column. 

Thus  the  numerator  is  kept  from  increasing  too  rapidly  so  that  inter¬ 
mediate  values  of  P  are,  except  in  unusual  instances,  kept  to  a 
value  less  than  one.  In  this  way  there  is  little  likelihood  of  exceed¬ 
ing  the  upper  limits  on  any  computer.  Never-the-less,  intermediate 
values  are  periodically  looked  at  and  adjusted,  if  necessary,  by 
the  scaling  procedure  described  below. 

ak  represents  each  cell  value.  It  is  set  at  one  and  is  increment¬ 
ed  by  one  up  to  the  value  of  a^  j,  reset  to  one,  incremented  up  to 
the  value  of  a  ,  etc.  Any  A.,  equal  to  zero  is  ignored  since  it  does 

lu  1J 

not  affect  any  marginal  totals  or  the  sample  size  and  its  factorial 

value  of  one  does  not  alter  the  calculations. 

d  represents  the  total  sample  size  and  is  incremented  by  one 

for  each  successive  calculation.  The  computation  is  complete  when 

a  has  been  incremented  to  the  final  value  of  r  . 

m 

The  final  value  of  P  must  be  less  than  one.  However,  it  is  this 
very  point  which  creates  the  most  formidable  problem  of  all.  As 
the  number  of  cells  and/or  the  sizes  of  the  marginal  totals  increase, 
the  number  of  points  in  the  sample  space  increase  astronomically 
while  corresponding  probabilities  get  unbelievably  small.  Con¬ 
sider  the  following  examples: 
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Table  IV 
An  Example 


2 

0 

3 

5 

±_ 

6 

5 

12 

3 

6 

8 

17 

(From  Ref  10) 


a.  The  probability  of  the  given  array 
is  .027149.  The  number  of  possible 
arrays,  given  the  marginal  totals, 
is  18. 


Table  V 

An  Example _  b.  The  probability  of  the  given  array 


2 

0 

1 

3 

is  .  023139. 

Although  the  sample  size 

1 

2 

2 

5 

2 

3 

4 

_2 

remains  the 

same,  the  number  of  pos. 

5 

5 

7 

17 

sible  arrays,  given  the  marginal  totals. 


is  179. 


Table  VI 

An  Example _  c.  The  probability  of  the  given  array 


11 

74 

181 

22 

288 

1 

25 

201 

109 

336 

12 

99 

382 

131 

624 

is  .77316  x  10  Among  the  171,  600 
possible  arrays  is  one  with  a  probability 


(From  Ref  5) 


of  .49981  x  10 


-94 


Quite  obviously  the  situation  exists  where  many  machines  would 
have  an  exponent  underflow.  Additionally,  in  trying  to  add  numbers 
such  as  these  a  sub-total  could  conceivably  be  reached  to  which 
each  additional  value  would  be  too  small  to  be  significant  and  would 
thus  be  lost  although  the  true  total  of  thousands  of  these  minute 
values  would  total  up  to  a  significant  sum.  The  problem  can  be 
overcome  by  scaling  each  value  and  considering  the  individual  prob¬ 
abilities  as  two-part  numbers.  Hence,  PSTOR(l)  could  represent 
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a  value  between  one  and  10  PSTOR(2)  would  represent  values 
-10  -20 

from  10  to  10  ,  and  PSTOR(S)  would  be  a  value  between 

lO’10^-1)  and  10  X  The  computer  merely  takes  a  value 

-47  40  -7 

like  10  ,  multiplies  it  by  10  ,  and  stores  it  as  PSTOR(5)  £  10 

Additions  to  PSTOR{5)  would  be  continued  until  the  total  exceeds 

-10 

one.  The  value  is  then  multiplied  by  10  and  put  in  location 
PSTOR(4).  At  the  completion  of  a  run,  the  lowest  ten  locations 
are  accumulated  and  PSTOR(l)  is  output  as  a  final  result.  Any 
values  remaining  in  other  locations  could  only  affect  the  result  by 
a  trivial  amount  and  can  be  ignored.  The  program  could  easily 
be  modified  to  change  scaling  factors  so  as  to  decrease  roundoff 
error  depending  on  the  machine  to  be  used. 

One  additional  step  is  required  when  this  first  probability  is 
computed.  It  must  be  used  for  comparison  with  all  other  probabil¬ 
ities  computed  later  on.  Some  decision  must  be  made  as  to  the 
amount  of  accuracy  desired,  consistent  with  machine  capabilities. 
The  present  program  was  tested  on  a  CDC  6600  with  29  significant 
digits  in  double  precision.  By  taking  the  initial  probability,  mul¬ 
tiplying  by  10*0,  truncating  to  an  integer  value,  all  subsequent 
probabilities  can  be  compared  to  10  decimal  places  of  accuracy. 

Any  probability  with  the  same  scaling  factor  as  that  originally  com¬ 
puted  is  similarly  manipulated  for  comparison.  This  procedure 
was  tested  with  example  b  above,  and  had  zero  roundoff  error. 


The  test  printed  out  all  possible  arrays  with  their  associated  prob- 


19 


GSA/MA/72-6 


abilities,  and  the  original  array,  reencountered  after  108  trials, 
showed  the  identical  probability  although  several  different  algorithms 
had  been  used  in  the  interim,  and  some  roundoff  error  would  have 
been  reasonable. 


Finding  All  Feasible  Arrays 

Within  the  constraints  of  fixed  marginal  totals  each  feasible 
array  must  be  determined.  Since  the  marginal  totals  are  fixed 
the  last  cell  in  each  row  can  be  determined  from  the  first  n-1 
cells  in  that  row  and  the  row  total.  Similarly,  the  last  cell  in  each 
column  can  be  determined  from  the  first  m-1  cells  and  the  column 
total.  Thus: 


and, 


a. 

in 


a. . 


for  all  i 


for  all  j 


(13) 


(14) 


Hence,  (n-1)  x  (m-1)  cells  and  the  marginal  totals  must  be  speci¬ 
fied  to  completely  determine  any  array.  Obviously  no  cell  value 
can  be  less  than  zero  nor  larger  than  the  smaller  total  for  the 
row  or  column  associated  with  that  cell.  The  bounds  on  any  cell 
value  are,  therefore, 


OS  aij  —  min(r  .,  Cj),  (i  =  1, 2 . m;  j  =  1,  2,  ....  n)  (15) 

Notice  that  the  (mn)  cell  is  somewhat  unique.  Its  value  can  be 
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0  to  Vj.  This  would  be  repeated  each  time  the  second  value  was 
increased  until  a  further  increase  would  exceed  v^.  At  this  time, 
could  be  set  to  one  while  and  were  set  to  zero  and  all  of 
the  preceeding  steps  would  be  repeated  for  each  increment  of  the 
third  value  until  it  could  no  longer  be  increased  without  exceeding 
its  limit  of  v^.  Thus,  (Vj  +  1)  x  (v2  +  1)  x  (v3  +  1)  possible  com¬ 
binations  would  be  generated.  This  is  essentially  the  procedure 
used  by  the  program  except  that  each  change  in  a  changeable  value 

must  be  followed  by  checks  to  determine  if  all  of  the  a.  's,  the 
'  in 

a  .'s,  and  the  a  are  positive, 
mj  mn 

Initially,  all  of  the  changeable  values  are  set  to  zero.  By 

Eq  (13)  a.  s  r.  (i  =  1,  2 . m-1),  by  Eq  (14)  a  =  c.  (j  -  1, 2,  . . . ,  n-1), 

ln  1  mj  J 

and  a-mn  is  then  computed  by  Eq  (17).  If  amn  is  non-negative  the 
array  is  feasible  and  its  probability  is  computed.  Whether  or  not 
amn  is  non-negative,  the  first  changeable  variable  (ajj)  *8  incre¬ 
mented  by  one,  a,  and  a  ,  are  each  decreased  by  one  and  a 
is  increased  by  one.  Any  time  these  last  three  values  are  non¬ 
negative  the  array  is  feasible.  Once  a  or  a  becomes  negative 
any  further  increase  in  a^  ^  is  futile  and  it  is  reset  to  zero.  The 
next  changeable  value  is  now  incremented  by  one  and  checks  for 
non-negativity  are  conducted,  and  ajj  >s  again  cycled  until  a^ 

or  a  goes  negative.  The  next  changeable  value  is  incremented 
ml 

again  and  the  necessary  checks  are  made.  Each  changeable  value 


in  turn  is  incremented  until  the  last  cell  in  the  associated  row  or 
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column  goes  negative.  The  incremented  value  which  caused  this 

is  then  set  to  zero  and  the  next  changeable  value  is  incremented 

by  one.  Notice  that,  except  for  a^j,  a  value  is  never  incremented 

unless  the  preceeding  value  is  set  to  zero  and,  also,  whenever 

a  value  is  incremented,  all  previous  values  must  be  at  zero. 

The  procedure  is  complete  when  the  last  changeable  value  cannot 

be  incremented  further  due  to  the  restriction  of  Eq  (15). 

The  program  takes  advantage  of  one  additional  relationship  to 

shorten  processing  time.  Any  time  a  single  a^  (i  -  1,2,...,  m- 1 ; 

j  =  1,  2, . . . ,  n-1)  is  incremented  the  only  other  changes  will  be  a 

decrease  in  the  value  of  the  associated  a.  and  a  .  and  an  increase 

>n  mj 

in  a  .  Therefore,  if  a  feasible  array  has  been  found  and  the 
mn 

next  increment  in  a  changeable  value  does  not  cause  a  negative 

value  for  the  associated  a.  or  a  .  the  new  array  is  feasible  and 

in  mj 

a  need  not  be  checked  since  it  must  still  be  non-negative.  (Con- 
mn 

versely,  a  can  never  decrease  in  value  unless  a  changeable  value 
mn 

is  reset  to  zero. )  The  complete  sequencing  for  the  example  from 
Table  IV  is  demonstrated  in  Table  VII  on  the  following  page. 

The  Probability  of  Each  Succeeding  Array 

The  first  feasible  array  encountered  has  its  probability  com¬ 
puted  using  the  same  procedure  developed  for  determining  the  prob¬ 
ability  of  the  given  array.  Each  successive  feasible  array  could 
use  the  same  algorithm  but  this  would  be  extremely  time  consuming 
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Table  Vll 

Sequencing  of  the  Example  From  Table  IV 


0) 

0 

0 

5 

(2) 

1 

0 

4 

(3) 

2 

0 

3 

(4) 

3 

0 

2 

3 

6 

3 

2 

6 

4 

1 

6 

5 

0 

6 

6 

(5) 

0 

1 

4 

(6) 

1 

1 

3 

(7) 

2 

1 

2 

(8) 

3 

1 

1 

3 

5 

4 

2 

5 

5 

1 

5 

6 

0 

5 

7 

(9) 

0 

2 

3 

(10) 

1 

2 

2 

(ID 

2 

2 

1 

(12) 

3 

2 

0 

3 

4 

5 

2 

4 

6 

1 

4 

7 

0 

4 

8 

(13) 

0 

3 

2 

(14) 

1 

3 

1 

(15) 

2 

3 

0 

3 

3 

6 

2 

3 

7 

1 

3 

8 

(16) 

0 

4 

1 

(17) 

1 

4 

0 

(18) 

0 

5 

0 

3 

2 

7 

2 

2 

8 

3 

1 

8 

(From  Ref  10) 


and  much  easier  methods  are  available. 

If  an  a_  is  incremented  by  one  and  the  array  is  feasible  then 

the  only  other  changes  in  the  array'  will  be  that  *mn  is  increased 

by  one  and  a.  and  a  .  are  decreased  by  one.  The  probability 
in  mj 

Pt+j  of  this  new  array  is  simply  the  previous  probability  Pt  multi¬ 
plied  by: 


(a  +l)(a  >1) 
in  mj 

(*;;>(«) 

mn 


(18) 


where  the  a's  are  current  values  in  the  new  array.  If  Pt  is  treated 

as  though  it  had  been  computed  by  Eq  (6)  then  Eq  (18)  reflects  the 

change  in  going  to  Pt  +  j  caused  by  the  change  in  four  cells  in  the 

denominator,  a41  and  a_  are  increased  while  a.  and  a  .  are 
>J  mn  in  mj 
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decreased.  This  decrease  is  taken  care  of  in  the  numerator  of  (18). 

If  more  than  one  of  the  a^'s  changes  a  similar  procedure  is 
used  to  compute  the  probability  of  the  next  feasible  array  starting 
with  the  previously  computed  probability.  Consider  first,  the  effect 


of  a  single  cell  which  can  decrease  in  value.  Let  a  •  s  a  -  d, 

pq  pq 

where:  a  *  is  the  new  cell  value,  a  is  the  old  cell  value  and  d 

pq  pq 

is  the  difference  in  their  values.  Then  1/a  !  would  be  equal  to 

pq 

l/(a  x  (a  -l)x‘*-x(a  -d+1)  x  a  *!)andl/a  *!  would  equal 

pq  pq  pq  pq  pq 

TT  v_o(apq_v)  *  (1/*  !)•  Hence  the  previous  probability  must  be 
TT  d  - 1 

multiplied  by  >1  y-o^apq'v^  to  account  for  the  effect  of  the  decrease 

in  a^.  If,  for  example,  a  probability  had  been  computed  with 

a__  =  15  and  in  the  next  feasible  array,  a  were  12,  then  the  change 
PS  pq 

for  this  cell  would  change  the  probability  by  a  factor  of,  15x14x13. 
Had  a  increased  by  d  to  a  value  a  *  then  the  old  probability  would 

pq  PS 

have  to  be  divided  by  TT^_j  (a^+v).  Thus,  if  more  than  one  of  the 
ajj's  (i  =  l,2,...,m-l;j  =  1, 2, . . . ,  n-1)  changes  the  new  and  old 
values  are  stored  as  well  as  the  new  and  old  values  for  any  changes 
to  any  a^  or  a^.  The  probability  for  the  previous  array  is  then 
adjusted  for  each  cell  which  has  changed  by  using  the  procedure 
just  discussed. 


The  probability  associated  with  each  feasible  array  is  now  com* 
pared  with  the  probability  of  the  given  array.  If  the  present  scaling 
factor  is  smaller  than  the  original  one,  the  associated  probability 
is  greater  than  that  of  the  given  array.  These  values  are  nover-the- 
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l«ss  accumulated  for  a  final  check  to  aee  that  all  probabilities  have 
been  determined  by  looking  for  an  overall  sum  of  one.  For  scaling 
factors  equal  to  the  factor  associated  with  the  given  array,  the 
digits  are  compared  as  previously  described.  Scaling  factors 
greater  than  that  associated  with  the  given  array  indicate  immedi¬ 
ately  a  lesser  probability. 


A  Partial  Numerical  Example 

Consider  example  b,  from  before  with  the  array  fron  Tr.ble  V 


2  0  1 

3 

1  2  2 

5 

2  3  4 

__9 

5  5  7 

17 

With  four  independent  variables  all  cell  values  can  be  determined 

from  the  marginal  totals  and  the  values  a  ,  a  ,  a  ,  and  a,.. 

11  12  21  “ 

The  probability  of  the  given  array  is  .023139  and  the  scaling  factor 
0 

is  10  .  a.  a ...  a  ,  and  a7_  are  set  to  zero  and  the  array  becomes 
11  12  21 


0  0  3 

3 

0  0  5 

5 

5  5-1 

9 

5  5  7 

17 

The  negative  entry  is  non-feasiblc  and  a^  is  incremented  by  one 
giving: 


1  0  2 

3 

0  0  5 

5 

4  5  0 

_2 

5  5  7 

17 

The  probability  for  this  array  is  .000077129.  This  is  less  than 
.023139  and  is  stored  in  PSTOR(l). 
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...  •  .  WI’  VIKW 


•  i«  inci 


sd  by  one  to  give: 


2  0  1 

3 

0  0  5 

5 

3  5  1 

9 

5  5  7 

17 

The  probability  for  this  array  is: 

P  s .000077129  x (2/2)  x  (4/1)  = .00030851 
This  probability  is  again  less  than  .023139  and  is  added  in  to  the 


value  in  PSTOR(l).  a^  is  again  increased  by  one  to  give: 


3  0  0 
0  0  5 
2  5  2 


3 

5 

9 


17 


5  5  7 

The  probability  for  this  array  is: 

P=  .00030851  x  (1/3)  x  (3/2)  =  .00015426 
After  computing  and  comparing  probabilities  the  next  increase  in 


j  gives: 


4 

0 

-1 

3 

0 

0 

5 

5 

1 

5 

3 

? 

5 

5 

7 

17 

The  negative  entry  is  non-feasible  so  a  ^  is  reset  to  zero  and  a^ 


is  incremented  by  one.  The  new  array  becomes: 


0  1  2 

3 

0  0  5 

5 

5  4  0 

9 

5  5  7 

17 

1 


Now  a  j  j  has  decreased  from  3  to  0,  a^  has  decreased  from  5  to  4, 


and  a  from  2  to  0.  a  increased  from  0  to  1,  a  from  0  to  2, 
3  3  1  13 
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and  a^  from  2  to  5.  The  new  probability  is: 

P=  .00015426  x  (3x2xl/l)  x  (5/2x1)  x  (2xl/5x4x3)  r  .000077129 
To  illustrate  other  intermediate  steps,  sequencing  for  feasible 
arrays  goes  from: 


0  3  0 

0  0  5 

'22 


0  0  3 
1  0  4 
4  5  0 


1  0  2 
1  0  4 

3  5  1 


and  from: 


0  3  0 
4  1  0 
1  1  7 


0  0  3 
0  2  3 
5  3  1 


1  0  2 
0  2  3 

4  3  2 


The  final  computation  is  on  the  array, 

3  0  0  3 

0  5  0  5 

2  0  7  9 

5  5  7  17 

The  final  outcome  showed  179  feasible  arrays.  Computation  was 
done  in  double  precision  with  results  printed  out  to  18  significant 
digits.  On  the  above  example,  the  probability  of  an  array  equal 
to  or  less  than  the  given  array  was  .776840806252570958.  The  prob¬ 
ability  of  an  array  more  probable  than  the  given  array  was  reported 
as  .223159193747429042.  Their  sum  was  1.000000000000000000. 


Final  Considerations 

As  alluded  to  previously,  the  number  of  feasible  arrays,  given 
fixed  marginal  totals,  increases  explosively  as  the  dimensions  and 


sample  size  of  any  given  table  increase.  The  number  of  possible 
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arrays  can  be  found  by: 


n1.n2.^..nk 

r-ik 

where  the  x-  are  cell  entries  in  the  range  of  0  to  n  and  )  .  n.  5  n 
1  L-ii-l  * 

(Ref  25:26,  27).  The  summation  in  Eq  (19)  is  then  taken  over  all 

sets  of  values  of  n^,  n^ . n^  such  that  their  sum  is  n  and  each  m 

is  an  integer  in  the  range  of  0  to  n  inclusive.  Finding  this  value  is 
a  formidable  task  in  itself  and  its  solution  does  not  provide  any 
real  guidance  in  estimating  computer  run  time.  Some  tables  se¬ 
quence  quite  readily  from  one  feasible  array  to  the  next  and  run 
time  is  relatively  short.  Other  tables  may  have  fewer  possibil¬ 
ities  but  require  much  more  checking  and  rearranging  with  a  re¬ 
sultant  increase  in  run  time.  As  a  consequence,  every  effort  was 
made  to  minimize  run  time  while  relying  on  the  FORTRAN  IV  lan¬ 
guage  (Ref  22)  which  is  applicable  to  most  machines.  As  a  result, 
time  saving  steps  available  to  FORTRAN  EXTENDED  (such  as 
mixed  mode  arithmetic)  were  purposely  avoided. 

The  program  w,  checked  on  a  CDC  6600  and  use  on  a  smaller 
machine  may  require  more  frequent  checks  for  scaling  than  is  now 
accomplished.  Additionally,  increased  accuracy  could  be  sought 
by  resorting  to  double  precision.  The  present  program  is  in  single 
precision  since  this  provided  a  savings  in  time  of  approximately  40%. 

The  CDC  6600  can  operate  in  the  range  of  10*^®  to  10 
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Real  constants  have  a  maximum  of  15  decimal  digits  while  double 
precision  constants  may  be  a  maximum  of  29  decimal  digits. 
Consequently,  use  of  the  program  on  a  smaller  machine  may  result 
in  exponent  underflows  or  loss  of  significance.  The  present  pro¬ 
gram  can  be  adapted  to  smaller  machines  with  minor  changes. 

A  Summary  Example 

Return  again  to  the  data  on  multiple  ejections  and  the  research¬ 
er's  results  following  Table  III.  The  computer  program  reported: 

1.  456  possible  arrays. 

\  V 

2.  The  probability  of  an  array  equal  to  or  less  than  the  given 
array  is  .  104957718232142. 

3.  The  probability  of  an  array  greater  than  the  given  array 
is  .895042281766887. 

4.  As  a  check,  the  sum  of  all  probabilities  is  .999999999999030. 
Hence,  a  "two-tailed"  test  of  independence  at  the  .  10  level  would 
fail  to  reject  and  no  statement  is  appropriate  {at  the  .  10  level) 
concerning  a  greater  or  lesser  likelihood  of  fatality  on  a  second 
ejection  following  non-injury  on  the  first  ejection. 
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V.  The  Resulting  Program 

Appendix  A  is  a  flow  chart  for  the  program.  The  program 
itself  is  included  as  Appendix  B.  Comments  are  included  so  that 
the  program  is  an  entity  in  itself  and  can  be  used  independently 
of  the  report  contained  here. 

Program  run  time  can  become  quite  large  and  is  dependent 
upon  the  number  of  rows  and  columns  as  well  as  the  sample  size. 
For  instance,  example  a  in  the  text  took  approximately  .2  seconds, 
while  example  b  took  .25  seconds.  Example  c,  with  a  fairly  large 
sample  size,  took  approximately  12.  seconds.  A  3  x  5  array  with 
all  row  totals  equal  to  five  and  all  column  totals  equal  to  three 
takes  about  one  second. 

A  final  note  of  caution  is  in  order.  The  table  must  be  arranged 
as  specified  in  the  comments.  Failure  to  do  so  may  return  incor¬ 


rect  answers. 
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VITA 


Marvin  Francis  Schwartz,  Jr.  was  born  I 


»> 


t' 


[Pll  Redacted] 


He  enlisted  in  the  United  States  Air  Force 

in  August  1950.  He  received  a  competitive  appointment  to  the 
United  States  Military  Academy  at  West  Point,  N.  Y.  and  entered 
with  the  Class  of  1956.  He  graduated  5th  in  a  class  of  480  and  was 
commissioned  in  the  United  States  Air  Force.  He  received  his 
pilot's  wings  in  August  1957  and  was  assigned  to  the  11th  Troop 
Carrier  Squadron  at  Dreux  Air  Base,  France.  Subsequently  he 
had  assignments  as:  Instructor,  Ground  Electronics  Officer;  and, 
Aide-de-Camp.  He  flew  593  combat  missions  as  a  Forward  Air 
Controller  in  Vietnam  and  followed  this  with  duty  as  a  KC-135  Air¬ 
craft  Commander.  Professional  military  education  includes 
Squadron  Officer  School  and  Air  Command  and  Staff  College.  He 
is  presently  assigned  in  the  Graduate  Systems  Analysis  course  at 
the  resident  school  of  the  Air  Force  Institute  of  Technology. 


Permanent  address: 


[Pll  Redacted] 
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